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Abstract 

We propose the generalized quadrature methods for numerical solution of singular integral equation of Abel type. We overcome the 
singularity using the analytic computation of the singular integral. The problem of solution of singular integral equation is reduced 
to nonsingular system of linear algebraic equations without shift meshes techniques employment. We also propose generalized 
quadrature method for solution of Abel equation using the singular integral. Relaxed errors bounds are derived. In order to 
improve the accuracy we use Tikhonov regularization method. We demonstrate the efficiency of proposed techniques on infrared 
tomography problem. Numerical experiments show that it make sense to apply regularization in case of highly noisy (about 10%) 
sources only. That is due to the fact that singular integral equations enjoy selfregularization property. 

Keywords: integral equations, singular kernels, quadrature, regularization, Abel equation, infrared tomography, midpoint 
quadrature. 


1. Introduction 


Numerical methods for solving a variety of singular integral equations (SIE) are offered in many publications, 
here readers may refer to l|2l-ll5l. llT?llT^IT8ll2Tll22ll^lZ71f3Tlf34lf37ll40ll4TI and others. A one-dimensional SIE of 
the 1st and and 2nd kind with Cauchy kernels Hilbert kernels, logarithmic et al., as well as two-dimensional, nonlinear 
SIE have been addressed. In present article we concentrate on Abel singular integral equation HI Em 111 1121111121 
|2l|26l[35l|36l|33[39ll40l 

R 

2 r k(r) dr - q(x). Q < x < R, (1) 

J Vr2 - 

A 

where k(r) is desired function, q(x) is the source function. Equations Q are widely used in practical models includ¬ 
ing plasma diagnostics, thermal tomography. X-ray CT, spectroscopy, galaxy clusters astrophysics, etc. In all these 
problems the object of interest enjoy the axial (or spherical) symmetry. Abel equation is also can be written as 


/ 


Hr) 

yjx- r 


dr — q(x), 0 < X < R. 


( 2 ) 


It has been studied in this form is Ill|5l[ni[ll[l8]|2ll|2ll2l|29l[3l]|3l|39llll]. Equation Q describes various 
problems in mechanics (such as tautohron problem), scattering and other problems. Of course, one may transfer SIE 
Q into SIE Q and vice versa, but it makes it more complicated to analyse their physical meaning. 

1 













V S. Sizikov, D. N. Sidorov / APNUM 00 (2015) /-fTT] 


2 


Let us below outline the main algorithms for numerical solution of SIE and singular integrals computation. 

For more details readers may refer to EglIS]. 

1. Algorithms based on relevant mesh shift. In papers Ejia the discrete meshes of knots with respect to variables 
r and x are introduced, i.e. rj - jh, Xj - ri + A, j,i - 0,1,, n, r„ - R, where step h =/?/«, A is mesh shift 
which is h/2 IlorAe (0,/z/2) H). Introduction of the shift A enables singularity overcome when it comes to 
quadrature rules application. But such algorithms need this shift selection. 

2. Quadrature type methods. One of the popular methods (here readers may refer to work |[3l) is Discrete Vortices 
Method where the integral with Cauchy kernel 


1 

In 


1 

/ 


y{x) 

X - Xq 


dx - f{xo), -1 < Xq < 1, 


is approximated with lift rectangles quadrature rule and using meshes on x and xq with shift A = hjl. This gives 
the system of linear algebraic equations (SLAB) with non zero main diagonal. 

In work fS) the “onion peeling” method for solution of SIE Q is suggested. Here region r 6 [0, R] is approxi¬ 
mated with rings Ar wide of constant values k e {rj - Ar/2, rj + Ar/2) for each rj. Here meshes are assumed to 

be uniform (A = 0). The main idea in this method is that integral y [ 2 C “ ^rjl > x, is computed 

analytically and its finite. Further midpoint quadrature is used resulting systems of linear algebraic equations 
with upper triangular matrix with respect to kj = k{xj). The similar method is suggested in lEl. 

3. Solution approximation. In works Il2ll24ll30l[3n[34ll40ll et al., the desired solution k(r) (as well as the right-hand 
side qix)) is approximated with an orthogonal polynomial, shifted Legendre polynomials, normalized Bernstein 
polynomials, algebraic or trigonometric polynomial or polynomial spline with coefficients determined with 
minimum of discrepancy between the left-hand side and right-hand side of Q. This leads to a projection 
method (the Galerkin method, the collocation method, the method of splines, the quadrature method, the least 
squares method, etc.) and to the solution of a SLAB wrt the corresponding polynomial coefficients. 

In these algorithms, there is a self-regularization, and in the case of using the relative shift of meshes, the shift 
A plays the role of the regularization parameter. Namely if A is closer to h/2 then solution k{r) is more stable, 
but it makes reduction of resolving capability of the method. If A is closer to zero, then solution is less stable 
but resolving capability of the method is higher. In all these algorithms, a SLAB is with prevailing (but not 
infinite) matrix diagonal. 

We also note a number of algorithms. Equation Q, as is known, has an analytical solution Ifni2ll6ll^ l24lf^[371 

R 

k(r) - — r ^ dx, 0 < r < R. (3) 

n J yj yf _ f.! 

r 

However, solution ([^ contains derivative q'{x) of experimental (noisy) function q{x) and the problem of dif¬ 
ferentiation is ill-posed 1381. Moreover, integral in Q is improper (singular). Nevertheless, a number of the 
following algorithms is proposed to compute the solution according to Q. 

4. Interpolation and quadrature method. In Q, derivative q'{x) was computed using interpolation on three (and 

two) neighboring points (discrete values of x). Integral £ ££2 computed analytically (without 

singularity). The similar algorithm was suggested in El using generalized left rectangles formula. 

5. Approximation of the right-hand side qix) is used in works HSl |24l ED). Function qix) is suggested to be 
approximated by a linear combination of smoothing polynomials (or splines) uniform for the whole interval 
X e [0, R]. Derivative q'ix) is computed using polynomial (or spline) differentiation. Solution k(r) in accordance 
with *0 is computed by summing the integral in Q along segments that performed analytically (see ll40l pp. 
188-189]). 

6 . Algorithm without using derivative q'ix). In 13 , formula Q is converted (by means of integration by parts) into 
the following expression that does not contain derivative q'ix) (cf. imsoi): 


kir) = - 


1 

n 


{q(R)-q{r) 1 

C X [qix) - qir)] ^ J 

1 — r^ 

j sjix^- r^f j 


2 


0<r<R. 
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This algorithm is implemented, e.g., in paper HOl pp. 217-220] using the cubic spline (see ll^ p. 273]) for 
q(x). 

7. Use of regularization. Abel’s equation 0 enjoys self-regularizing property due to the singularity, as a result the 
problem of its solving is moderately ill-posed |8]. This means that above mentioned algorithms are moderately 
stable. Nevertheless, in papers lU 0 HI US et ah, the Tikhonov regularization method IfTOl fT6][38l ? ] was used 
to enhance the stability of algorithms. 

In this work, we develop the following variant for numerical solving some SIE. We make the meshes of 
nodes in r and x coincide (i.e. A = 0) and eliminate the singularities using the generalized quadrature formula (cf. 
ifTOilMllJTll l. However, such a technique can be applied to not all SIE. Eor example, it is not applicable to SIE with the 
Cauchy kernel, but it is applicable to some SIE with logarithmic and other (weakly singular) kernels. In this paper, we 
consider the solution of equation ([1]) wrt the desired function k(r), as well as numerical computation of k(r) according 
to ([^ by the generalized quadrature method with use of Tikhonov regularization. 

2. The generalized quadrature method 

Let us describe method using generalized left rectangles formula in application to numerical solving SIE Q 
(the first method) and to computation of k{r) according to ([^ (the second method). 

It is to be noted here that in ||6l|8l, the numerical method “onion-peeling” is suggested for computation of 
integrals in ([1]) |[8l and in ([^ jh]. Here, the uniform coinciding node meshes in r and x and the middle rectangles 
quadrature formula have been employed. In llTTlI . also uniform coinciding meshes have been employed combined 
with more usable the left rectangles formula. 

In present paper, we use nonuniform meshes and left rectangles resulting more generic and convenient al¬ 
gorithm. The solution error estimates for equation ([1]) by the generalized quadrature method are also derived. This 
method is described below in two variants (the first and second methods). 

2.7. First quadrature method 

Eirst quadrature method employs generalized left rectangle formula. Let us introduce nonuniform (but coin¬ 
ciding) meshes on x and r as follows 

Q - X\ - r\ < X2 - r2 < ■ ■ ■ < Xi — ri < ... < Xn — r„ — R. ( 4 ) 

Here, R - rmax is boundary value such as k{R H- 0) = 0. On each interval {rj, r^+i), j - 1,2,...,«- 1 we suppose 
approximately 

k{r) - k(rj) = kj - const. (5) 

We have the following 

Lemma 2.1. Under condition (|^, one has the equality 

0+1 _ 

n 

ye[l,n-l], x < 



Proof. Integral (table) 


/ 


: dr - — x^ for x < r. 


□ 


VjTT 

whence, taking into account (5), we obtain (6). 

Definition 2.1. We call formula 0 as generalized quadrature formula of left rectangles (cf 1 1791? ) for the specific 
singularity rj — x^, and multipliers — x^ — ^r^. — x^ are the quadrature coefficients of this singularity. 
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The specifics of the formula (|^ is that the singular integral £ '*' dr is calculated analytically accurate 

and without peculiarity. If it is calculated numerically by the usual left rectangles quadrature formula, then at x - rj 
there will be a division by zero. 

Let us now formulate the main result as following 


Theorem 2.1. Numerical solution of equation Q according to the first quadrature method is defined as following 
recursion 

Q^n-l 


^n-\ 


Pn-\,n-\ ’ 


kj = 


qi 


Z n-i 


Pij 


i — n — 2,n — 3,... 


(7) 


where ki = k{ri), qi = qixj). 

Proof. Integral in Q is sum of integrals (|^, i.e. 


kn = ^„-2 + { r[\ J’PU (kn-1 - kn^2), 

pu = - 4 


r^. - x^, 

j I ’ 


(8) 


i\ ^ j 


(9) 


This is the SLAB wrt {kj}"^^. SLAB (0 is upper triangular and its solution can be recursively constructed. Brom (0 
for ! = n - 1,n - 2,..., 1 we eventually obtain k„_i, k„^ 2 , ■ ■ ■ ,ki according to Q. As to the value of k„ = k{R), it can 
not be found by this scheme, but can be additionally determined as k,, - 0 from physical concepts or kn - kn-i or can 
be derived using linear extrapolation iTSl . as was done in 0- □ 

In Is], formulae of type (|^ are also given, but for the case of uniform (and coinciding) meshes in r and x and 
using the middle rectangles quadrature formula. Burthermore, important formulae of type 0 is not given. 

Bormulae Q-(|9]) are more common and more convenient than in |[8l and formulae (|7]) give a solution in the 
explicit form. The method according to Q-Q for solving the equation ([1]) is called in Il37l the generalized quadrature 
method for solving SIB Q. This paper presents a more general formulae Q and Q than in llJTl . In Sec. 3, estimates 
of the errors for this method are given. 


2.2. Second quadrature method 

The second method is generalized quadrature method for computation of singular integral 0 giving the 
solution k{r). Let us assume k'{x) to be computed with some stable method. As in the first method we introduce node 
meshes 0. On each [x;, x;+i), i = 1,2,..., n - 1 we assume 


q'(x) - q'(xi) = q'l - const. 


Let us fomulate the following 

Lemma 2.2. Under condition (|10[i the following equality is true 


Xj-H 

f 


q'(x) ■^'■+1 

■.dx = \n 


4 




Vx^ - 


, + 


i — 1,2,... ,n — I, r < Xi < Xi+\ < R. 


( 10 ) 


( 11 ) 


Proof. Integral J ^=3 = In (■* + Vx^ - ^ for r < x is the table integral. Taking into account the condition ( [TOl i we 

obtain ([n}. □ 


4 
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Definition 2.2. Formula (|1 l|i is generalized left rectangle quadrature rule for the singularity 1 / slr^ — x^, and multi¬ 


pliers In 




JT- 


are quadrature coefficients of this singularity. 


Now we can formulate 

Theorem 2.2. Numerical solution ofSIE Q by formula 0 according to the second generalized quadrature method 
is result of the following recurrence formulae 


where 








]> '^j 


n 


Xiv\ 


gij 


= In 


- 


Xi + 




Proof. Integral in ([^ 


IS sum 

R 


f 4 


of integrals o over the separate intervals [xi, Xi+i), i.e. 

n-2 

dx—'^^\\\ - ^ - q'-. j — 1,2,... ,n — 1. 


q'(x) 


x^-r^. 


Xi+l + 


X,-+ ^ 

Ixf-rj 


( 12 ) 


(13) 


(14) 


As a result, solution 0 in the discrete form is according to ( [T^ , ( [T3] l. For j = I this (second) method due to 

( [T4| ) gives uncertainty oo ■ 0 for i = y = 1 since xi - ri - q'^ - 0. In this case let’s use the hrst method to determine 
ki. Using 0 and ([^ for i — 1 we hnd ki, ref. ( |12| i. As to kn it can be calculated using linear extrapolation (see Q). 

□ 

The advantage of the above two methods is that they do not require the relative shift of meshes and their 
integrals with singularities rj - x^ and 1/ Vr^ - x^ are calculated analytically and without divergences. 

However, these methods are not suitable for all singularities, e.g., for numerical computation of hypersingular 
integral with the Cauchy kernel ^dr^^. 


3. Error estimates 

Let us derive errors estimates for solution of SIE ([T]) using hst quadrature method (cf. llJ7l[39ll20l '). 


3.1. Quadrature error on small interval 

Let us estimate quadrature error for computing integral (|^ on separate small interval {rj, rj+i) due to approx¬ 
imation 0 (while without the measurement error for q{x)). 

Lemma 3.1. Integral 

0+1 

r 

k(r)dr, Xi < r,- < r,+i < R, 

(15) 

/ = 1,2,...,n — 1, j — i,... ,n — I, 

when using the generalized left rectangles formula 0 and taking account of quadrature error caused by the approxi¬ 
mation 0 is equal to (refinement of formula 0) 


f 


I 



k(r) dr - pqkj -t- Asq, 


(16) 


5 
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where pij are quadrature coefficients and Asij is quadrature error of computation of integral •EH) approximately 
equal 




(17) 


fi 6 [rprj+f). 


fj 6 [O’O+i)- 


Proof. Using the first method we assume k(r) = kj, r e [rj,rj+i) (see (0), i.e. we represent function k{r) by the 
interpolation Lagrange zero degree polynomial lEl. Error of such interpolation is Akfr) = k(r) - kj = k'{f) (r - rf), 
where f - fj(r) e [rj,rj+i). Then 

k(r)^kj + k'(f)(r-rj). (18) 

Let us now substitute El into El, we get ([Thll, where 



(19) 


Integral in ( [T9] l can be analytically computed giving us an estimate ([T^. □ 

It is to be noted that derivative k'{ff) in ([TtIi can be approximated with 



( 20 ) 


or by other way ET). 

3.2. Quadrature error 

Let us estimate quadrature error of the solution of equation 0 due to the approach 0 (without error qix)). 
We formulate this as a theorem. 

Theorem 3.1. Errors of numerical solution of equation 0 by the first generalized quadrature method according to 
0 are computed with the following recurrence 


Ak„ — Ak„^i — 




( 21 ) 


where 


n-\ 



( 22 ) 


Here, pij, Asq and k'iff) are computed based on 0,E1 and respectively. 
Proof. Let us write integral in 0 as sum of integrals wrt intervals 




(kj + k\f)(r-rf})dr, i = 1,2,1. 


6 
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Then 


«-i 

5 / 


0 a/" 


: lS.kj{r) dr 


n I '■^+‘ 

_v fo) 




k'{^)dr, i — \,2,... ,n — 


(23) 


In order to compute the integral in the left-hand side of ( |2^ , we employ the formula of left rectangles, i.e. 
we assume Akj(r) - ISki/j) = Akj - const, r e [r^, r^+i) and compute the integral using the generalized formula in 
similar way with Integral in the right-hand side of ( |23| ) is equal to Asij due to ( |19| ). Then 

n-l 

^ Pij Akj — Si, i — 1,2,... ,n — 1, (24) 


7=' 


where e, denote sum ( |22| l. Here, ( |24l l is a SLAB wrt {Akj]"J^. It is also assumed that are computed using Q 
in advance. Then values are computed using ( |22] i, ([T^ and ( |20| ). We solve SLAB ( |24l i with upper triangular 

matrix and obtain solution ( |2T| l in the recurrent form, adding the condition Ak„ - Ak„^i. □ 

Remark 2. It is to be noted here that formulae ( |2T| l give errors of the solution Aki with their signs (cf. llJTl [39)) 
in contrast with other works where absolute values |AA:,| or upper bounds |AA:,| < ... or upper bounds by the norm 
||Ak,|| < ..., etc. are given. This enable us to obtain the refined solution 

ki - ki + Akj, i - 1,2,..., n, (25) 


using from Q and errors {Ak,)"^j from ( |2T| i. 

Remark 3. Brrors of numerical solution given in (|2TJ are obtained with regard to only the quadrature errors and 
the error of the right-hand side q(x) of equation Q is set equal to zero. Let us take into account the measurement 
errors of the source function qix). In Il20ll39ll . the error estimates for numerical solution of the Volterra integral 
equations of the first and second kind are derived taking into account both quadrature and source function errors. In 
similar way we can generalize recurrence formulae (|2T]l to the case of errors as follows 




|Ak„| — |Ak„_i| 

|e/l + 6i + Pij\Akj 


Pn-l,n-l 
^ n-l 


\Aki\ = 


i — n — 2,n — 3,... ,1. 


But estimates 


give overstated estimates of {|Ak,|)"^j due to using the operation | ■ | (absolute value). 


(26) 


4. Numerical illustration 


4.1. Software implementation 

Proposed two generalized quadrature methods have been implemented in MatLab 7.10 (R2010a). Bollowing 
the first method we search using Q, 0. Brrors {AkilJLj are calculated using ( |21[ ), as well as ([^, ( |17| i, ( |20| i and 
( |22| ). Refined solution is calculated with ( |25] l. Tikhonov regularization IfTOl[T6ll38l[^ 

k„ = {aE + A^AY^A^f (27) 

is employed, where the discrepancy principle ESll is used for choosing the regularization parameter a > 0: 

||Ak„-/|| = A (28) 

Here, / = ql2, E is identity matrix, A is matrix of the SLAB (0 represented as 

Ak = /, (29) 


where 


Aij = 


Pij, i ^ i, 

0, otherwise. 


/, 7 = 1,2,..., n — 1. 


(30) 


Bollowing the second method, solution k{r) has been computed using singular integral in (j^ based on gener¬ 
alized left rectangles formula according to ( [T^ and ( [T3| l. 

The first method was developed and implemented in software in more detail than the second method. 


7 
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4.2. Infrared tomography example 

The first method has been applied to axially symmetric flame diagnostics using infrared tomography iniisiiTi 
Elm nails] ESI Eg). Fig. shows measured output intensity I^nix) of rays (m from measurement) which go through 
gas, undergo absorption and emission and are accepted by detectors. Measurements have been performed in Technical 
University of Denmark, Department of Chemical and Biochemical Engineering (before 1 January 2012 Ris0 DTU) 
within a joint project IfTafTTI . 



Figure 1. Measured noisy intensity (difference of intensity in active and passive regimes). The mesh is nonuniform, number of nodes n = 11. 

Intensity I^ix) is recalculated into q„f x) - - ln[/ni(jii:)/B(7’o)] (the right-hand side of equation ([T]i), where 
B{Tq) is the Planck function of rays source with its temperature Tq - 894.4°C. Fig. j^shows function qaf x). 



Figure 2. Dimensionless right-hand side qmix) of SIE jl], « = 11. 

Fig. shows results of solution of SIE Q wrt absorption coefficient kj^ir) by the first method of generalized 
quadratures according to 0 and 0. 

As we see, the solution k^ir) suffer from significant artificial perturbations. This is due to too great step of 
node mesh (in other words, the smallness of n), as well as measurement errors in function Imix). Here we also demon¬ 
strate behavior of solution kma(r) derived with Tikhonov regularization using ((0, Q and 0. Regularization 
parameter a is chosen using discrepancy principle ( |28| ) where 6 - 0.037, as a result a — 10 = 0.813. Fig. 

demonstrates that solution has been smoothed by regularization method. 

To reduce the grid step in x as well as to moderately smooth the fluctuations in the function /m(x), a spline 
approximation was used lT5l l40l . Fig. 4 shows an approximation of the function /m(x) by cubic smoothing spline 
using the m-function csaps.m. 

The smoothed values of I{x) were generated with splines (Fig. 0 and then SIE 0 was resolved with gen¬ 
eralized quadrature 0. Fig. 0shows obtained solution k{r). We also applied solution using Tikhonov regularization 
( |27| ) for a - 10^^. Fig. 0 shows regularized solution ka{r). 

8 
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Figure 3. Absorption coefficient computed by the first generalized quadrature method and /:ma(^) computed by Tikhonov regularization, 



Figure 4. Measured Im{x) values are marked with o (n=l 1), values for n = 20 marked with • and spline interpolated values are marked with 
solid line. 



Figure 5. Absorption coefficient k{r) (without regularization) and ka{r) (with regularization) after spline smoothing of /m(-^), = 20, cm 


9 








V S. Sizikov, D. N. Sidorov /APNUM 00 (2015) j-fTT] 


10 


Fig. |^and|^demonstrate that application of spline based smoothing enable mesh step reduction for x (causing 
increase of n). This allows (moderate) smoothing k{r) and ka(r). Moreover, in case of noisy I{x) and big step of the 
mesh, regularization slightly improves solution as shown in Fig. In case of < 1% errors and small step (Fig. 
solutions k{r) (without regularization) and k„(r) (with regularization) are obtained practically the same (Fig|^. It 
confirms that the problem of solving the singular integral equations is moderately ill-posed and has the property of 
self-regularization. 


5. Conclusion 

In this paper we outlined two new methods of numerical solution of singular integral equation (SIE) of Abel 
type. The methods are based on the use of generalized quadrature formula of left rectangles. Specificity of methods is 
that singular integrals are computed analytically and without peculiarities. We derive recurrence formulae for solution, 
generally speaking, on a nonuniform node mesh. Estimates of quadrature errors of solution with regard to their sign 
in the absence and in the presence of errors in the right-hand side are found. In order to enhance the stability of the 
solution we used Tikhonov regularization. However, SIE enjoy self-regularization, therefore it is advisable to apply 
the Tikhonov regularization method only if there is a significant error (~ 10%) in the right-hand side and rough mesh 
step (when the number of nodes is small; n ~ 10). The method has been applied for solution of infrared tomography 
problem. 
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